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HIGHLIGHTS 


•  A  3-D  lattice  Boltzmann  model  of  cathode  GDL  in  PEM  fuel  cell  is  presented. 

•  Non-homogeneity  and  anisotropy  of  GDL,  and  electrochemical  reaction  are  included. 

•  General  orientation  of  GDL  carbon  fibers,  shapes  species  and  current  distributions. 

•  Normally  oriented  carbon  fibers  cause  larger  current  variations  on  catalyst  layer. 

•  Microstructure  must  be  regarded  as  an  important  characteristic  of  any  GDL  material. 
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High  power  density,  low  operation  temperature,  high  efficiency  and  low  emissions  have  granted  proton 
exchange  membrane  (PEM)  fuel  cells  the  most  promising  future  among  all  types  of  fuel  cells.  The  porous 
electrodes  of  PEM  fuel  cells  have  a  complicated,  non-homogeneous,  anisotropic  microstructure.  Therefore, 
pore-scale  modeling  techniques  such  as  lattice  Boltzmann  method,  which  can  capture  non-homogeneous 
and  anisotropic  microstructures,  have  recently  gained  a  great  attention.  In  the  present  study,  a  three- 
dimensional  lattice  Boltzmann  model  of  a  PEM  fuel  cell  cathode  electrode  is  proposed  in  which  electro¬ 
chemical  reaction  on  the  catalyst  layer  and  microstructure  of  GDL  are  taken  into  account.  The  model  enables 
us  to  simulate  single-phase,  multi-species  reactive  flow  in  a  heterogeneous,  anisotropic  gas  diffusion  layer 
through  an  active  approach.  To  show  the  capability  of  the  proposed  model,  reactive  flow  in  three  recon¬ 
structed  GDLs  with  different  anisotropic  characteristics  is  simulated  to  investigate  the  effects  of  GDL 
microstructure  on  species  and  current  density  distributions.  The  results  demonstrate  that  when  carbon 
fibers  are  more  likely  oriented  normal  to  the  catalyst  layer,  species  density  distribution  is  thicker  and  more 
disturbed.  Current  density  also  experiences  a  larger  variation  on  the  catalyst  layer  in  such  a  case. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  sustainable,  high  quality  life  is  the  basic  driver  for  providing  a 
clean,  safe,  reliable  and  secure  energy  supply  in  the  world.  Fuel  cells 
and  hydrogen  systems  are  important  technologies  which  may 
contribute  positively  to  the  future  world  energy  demand.  Among 
different  types  of  fuel  cells,  proton  exchange  membrane  (PEM)  fuel 
cells  currently  hold  the  most  promising  prospect.  They  can  be 
employed  in  a  wide  range  of  applications,  ranging  from  very  small 
portable  devices  such  as  mobile  phones  and  laptops,  through 
transport  applications  like  cars,  delivery  vehicles,  buses  and  ships, 
to  combined  heat  and  power  generators  in  stationary  applications 
in  the  domestic  and  industrial  sectors  [1]. 
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Nevertheless,  some  barriers  such  as  high  cost  and  operational 
difficulties  hinder  PEM  fuel  cells’  wide  commercialization,  which 
make  fundamental  research  for  their  development  inevitable. 
Modeling  of  reactant  gas  transport  accompanied  by  electrochemical 
reactions  in  the  electrodes  is  of  great  importance,  especially  in  the 
cathode  where  the  oxygen  reduction  reaction  is  sluggish  and  inef¬ 
ficient  [2].  Numerous  numerical  models  of  the  cathode  electrode 
with  different  features  can  be  found  in  the  literature  [3].  In  most  of 
these  models,  the  gas  diffusion  layer  (GDL)  of  cathode  electrode  is 
considered  as  a  homogeneous  and  isotropic  porous  medium,  while 
in  reality  that  is  not  the  case  [4].  In  this  regard,  pore-scale  modeling 
techniques  which  can  realistically  capture  the  microstructure  of  the 
cathode  GDL  are  extremely  attractive. 

Some  researchers  have  used  more  conventional  methods,  such  as 
the  finite  volume,  to  model  electrodes  of  a  PEMFC.  Harvey  et  al.  5],  for 
example,  used  such  a  method  to  investigate  different  approaches  to 
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Nomenclature 

w 

weighting  factor 

a 

roughness  factor 

Greek  symbols 

~c 

particle  velocity  (lu  ts-1) 

a 

transfer  coefficient 

cs 

speed  of  sound  in  lattice  (lu  ts-1) 

V 

activation  over-potential 

D 

diffusivity  (m2  s-1) 

A„ 

Eigen-parameter  in  Eq.  (11) 

dlb 

lattice  Boltzmann  diffusivity 

F 

dynamic  viscosity  (kg  m-1  s-1) 

Da 

Damkohler  number 

V 

kinematic  viscosity  (lu2  ts-1) 

F 

Faraday’s  constant  (A  s  mol-1) 

P 

density  (lm  lu-3  or  kg  m-3) 

Fo 

Fourier  number 

P* 

dimensionless  density 

f 

density  distribution  function 

T 

relaxation  time  (ts / 

j 

current  density  (A  m-2) 

k 

rate  constant  (m  s-1) 

Subscripts  and  superscripts 

kLB 

lattice  Boltzmann  rate  constant 

A 

type  A  species 

l 

characteristic  length  scale  (m) 

C 

type  C  species 

MW 

molar  mass  (kg  kmol-1) 

eq 

equilibrium 

P 

pressure  (Pa) 

i 

direction  i  of  lattice 

Ru 

universal  gas  constant  (J  mol-1  K-1) 

k 

species  k 

y 

particle  position  vector  (lu) 

LB 

lattice  Boltzmann 

r" 

reaction  rate  per  unit  surface  area  (mol  m-2  s-1) 

n 

nitrogen 

T 

temperature  (K) 

0 

oxygen 

t 

time  (ts) 

ref 

reference 

T t 

velocity  vector  (lu  ts-1) 

sr 

surface  reaction 

It 

composite  velocity  vector  in  Eq.  (5)  (lu  ts-1) 

w 

water 

modeling  the  catalyst  layer  of  PEM  fuel  cells.  They  compared  three 
catalyst  treatments  of  thin-film  model,  discrete-catalyst  volume 
model  and  agglomerate  model,  and  concluded  that  only  the 
agglomerate  model  could  capture  mass  transport  limitations  that 
occur  at  high  current  densities.  Sahraoui  et  al.  [6  presented  a  two- 
dimensional,  finite-volume  model  of  a  PEM  fuel  cell  taking  into 
consideration  the  finite  thicknesses  of  the  catalyst  layer  and  mem¬ 
brane.  Sun  et  al.  7]  also  introduced  a  two-dimensional  model  for  PEM 
fuel  cell  cathode  that  treated  the  catalyst  layer  as  agglomerates  of 
polymer  electrolyte  coated  catalyst  particles.  They  showed  that  the 
cathode  over-potential  inside  the  catalyst  layer  was  non-uniform 
influenced  by  the  channel-land  geometry. 

In  the  past  decade,  lattice  Boltzmann  method  has  emerged  as  a 
powerful  pore-scale  modeling  technique  applicable  in  complicated, 
heterogeneous,  anisotropic  porous  media  such  as  GDLs  of  PEM  fuel 
cells.  Lattice  Boltzmann  method  is  superior  to  other  conventional 
numerical  modeling  techniques  through  many  aspects  such  as 
capability  of  dealing  with  complex  boundaries  of  a  complicated 
microstructure,  easy  parallelizable  algorithm  development,  and 
facility  in  modeling  multi-phase  fluid  flow  in  a  porous  medium  [8]. 

However,  a  significant  challenge  for  lattice  Boltzmann  modeling  of 
the  cathode  electrode  is  a  proper  way  to  consider  electrochemical 
reaction  in  the  catalyst  layer  [9].  To  the  best  of  the  authors’  knowledge 
the  electrochemical  reaction  is  taken  into  account  in  only  a  few  two- 
dimensional  lattice  Boltzmann  investigations,  such  as  those  pre¬ 
sented  by  Chen  et  al.  [10—12];  moreover,  no  such  three-dimensional 
modeling,  taking  into  account  the  electrochemical  reaction,  was 
found  in  the  literature.  In  the  present  study,  a  single-phase,  three- 
dimensional  model  of  the  cathode  electrode  of  a  PEM  fuel  cell  is  pre¬ 
sented  in  which  the  catalyst  layer  is  considered  as  a  thin  interface 
where  the  electrochemical  reaction  occurs.  The  proposed  model  is 
implemented  to  investigate  the  effects  of  GDL  microstructure  on  the 
distributions  of  species  and  current  density.  An  advantage  of  the 
proposed  model  over  previous  lattice  Boltzmann  models  [10-12]  is 
the  use  of  an  active  approach  for  species  transport.  In  the  preceding 
investigations  [10-12],  a  passive  approach  was  applied  in  which  the 
velocity  field  is  only  solved  for  nitrogen  that  is  taken  as  the  solvent. 
This  may  lead  to  lack  of  accuracy  especially  when  the  density  of  other 


species  is  comparable  to  that  of  nitrogen  which  may  happen  at  high 
over-potentials;  in  such  cases  no  species  can  be  considered  as  a  solvent 
[13].  On  the  other  hand,  in  an  active  approach  the  internally  coupled 
velocity  fields  of  all  species  are  solved  individually  which  can  lead  to 
more  accurate  results. 

2.  Numerical  method 

2.1.  Framework  of  lattice  Boltzmann  method 

Lattice  Boltzmann  method  is  a  powerful  numerical  modeling 
technique  that  has  attracted  much  interest  in  the  past  decade.  In  this 
method,  a  fluid  is  assumed  to  be  composed  of  virtual  fluid  particles 
which  move  and  collide  with  each  other  in  a  lattice  structure.  These 
particles  are  treated  by  a  distribution  function  rather  than  by  posi¬ 
tions  and  velocity  vectors  [14].  This  kinetic  nature  of  lattice  Boltz¬ 
mann  method  provides  some  superior  features  relative  to  other 
numerical  methods,  such  as  the  linearity  of  the  transport  equation 
(in  comparison  with  nonlinear  Navier- Stokes  equations),  calcula¬ 
tion  of  pressure  by  an  equation  of  state  (in  comparison  with  calcu¬ 
lating  the  pressure  by  solving  Poisson  differential  equation),  etc.  [8]. 

The  present  numerical  model  is  based  on  the  lattice-Boltzmann 
method  with  a  single  relaxation  time  collision  operator  (the  so- 
called  BGK  model  [15])  and  the  popular  D3Q19  lattice  scheme. 
The  3  in  D3Q19  refers  to  the  number  of  lattice  dimensions  and  the 
19  denotes  the  number  of  possible  directions  for  a  particle  move¬ 
ment  in  the  lattice  as  shown  in  Fig.  1. 

The  lattice  Boltzmann  equation  is  derived  from  simplification  of 
Boltzmann  equation  in  a  lattice  [16].  It  can  be  expressed  as 

fiCr  +  ~CiM,t  +  M)  =fi(r,t)  +  Y[fleqC?,t)-fi(r,t)\  (1) 

where/]  is  the  density  distribution  function  in  direction  i,/-eq  is  the 
equilibrium  density  distribution  function,  T  refers  to  space  posi¬ 
tion,  t  is  time,  "cV  is  the  particle  velocity  vector  in  direction  i,  and  t 


1  lu  and  ts  are  the  unit  of  length  and  time  in  lattice  Boltzmann  method,  respectively. 
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is  the  relaxation  time  which  is  related  to  kinematic  viscosity,  v,  by 
t  =  3v  +  1/2.  The  equilibrium  density  distribution  function  is 
determined  by 


r 


=  WiP 


Cyt  1  (' CyU )2 

C}  +2  Cf 


IT t-Tt 

2  ~cf 


(2) 


In  this  equation,  w*  is  the  weighting  factor,  p  =  is  the  fluid 
density,  7?  =  is  fluid  velocity  and  cs  U  the  speed  of 

sound  in  the  lattice.  1 

Eq.  (1)  is  solved  through  collision  and  streaming  processes 
which  are  shown  for  D3Q19  scheme  in  Eqs.  (3)  and  (4)  respectively. 


fi(x,y,z,t  +  At)  =  fi(x,y,z,  t) 


1-M 

T 


+  v-C(w,  t) 


(3) 


/i(x  +  Ax,y  + Ay,z  + Az,t  + At)  =  fi(x,y,z,  t  +  At)  (4) 

To  completely  determine  density  distribution  functions,  f  s 
must  be  specified  at  the  boundaries.  A  simple  and  powerful  LBM 
boundary  condition  applied  to  no-slip  walls  is  called  bounce  back 
boundary  condition.  This  boundary  condition  is  based  on  the  idea 
that  particles  colliding  a  wall  in  a  direction  will  bounce  back  in  the 
opposite  direction  [17].  In  fact,  this  boundary  condition  enables 
LBM  to  model  fluid  flow  in  geometries  with  complicated  micro¬ 
structure  such  as  pore  space  of  a  porous  medium.  Several  versions 
of  this  boundary  condition  have  been  proposed  [13];  however,  the 
mostly  used  version  is  half-way  method  in  which  the  wall  is  located 
half-way  between  two  neighboring  grids. 


2.2.  Simulation  of  single-phase  multi-component  flow 

Generally,  two  passive  and  active  approaches  are  used  for 
simulating  single-phase,  multi-component  fluid  flows.  In  the  pas¬ 
sive  approach  the  chemical  species  whose  mass  fraction  is  greater 
than  the  others  is  considered  as  a  solvent  while  other  species  are 
considered  as  solutes.  In  this  way  the  velocity  field  is  only  solved  for 
the  solvent  species,  followed  by  the  solution  of  advection— diffusion 
equation  for  the  other  species.  Therefore,  the  velocity  field  is  only 
dependent  on  the  solvent  species.  When  species  concentrations  are 
comparable  and  no  species  can  be  considered  as  a  solvent,  this 
approach  may  lead  to  lack  of  accuracy  [13]. 

On  the  other  hand,  in  the  active  approach  the  internally  coupled 
velocity  fields  of  all  species  are  solved  individually;  hence,  more 
accurate  results  can  be  obtained  by  the  more  realistic  active 
approach.  To  apply  this  approach  in  lattice  Boltzmann  method, 
multi-phase  multi-component  models  are  adopted,  and  some  of 
their  settings  are  modified  to  reduce  the  number  of  phases  to  one. 

Several  multi-phase  multi-component  lattice  Boltzmann  models 
[18-20]  have  been  presented  in  the  literature.  In  the  present  study, 
the  Shan  and  Chen  model  [  19  ]  is  selected  over  other  models  due  to  its 


simplicity  and  easy  implementation.  Since  the  intention  here  is  to 
model  a  single-phase  reactive  flow,  the  number  of  phases  is  reduced 
to  one  simply  by  eliminating  the  inter-particle  forces  [  13  ].  Therefore, 
to  implement  Shan  and  Chen  model  for  the  single-phase  flow  con¬ 
sisting  of  oxygen,  nitrogen  and  water  vapor  species,  streaming  and 
collision  processes  are  solved  individually  for  each  species  to  find  fk, 
where  the  superscript  k  denotes  the  kth  chemical  species.  However, 
in  collision  process  instead  of  7 1  in  Eq.  (2),  the  composite  velocity, 
7 f,  must  be  incorporated  to  couple  the  species  lattice  Boltzmann 
equations.  The  composite  velocity  is  defined  as 


u  = 


E? 

k  i _ 

£>fc 

k 


(5) 


In  the  above  equation  pk  =  is  the  density  of  species  k  and 
xk  is  the  relaxation  time  of  species  k.  This  relaxation  time  is  related 
to  the  kinematic  viscosity  of  species  k,  vk ,  by  the  relation 

Tk  =  3vk  +  1/2. 


2.3.  Surface  reaction  model  based  on  modified  bounce  back 
boundary  condition 


Kamali  et  al.  [21]  recently  proposed  a  half-way  bounce  back 
boundary  condition  for  simulating  surface  reaction.  Their  proposed 
boundary  condition  is  only  applicable  to  a  typical  surface  reaction 
A  +  B  + ...  — ►  C  +  E  + ...  which  is  first  order  in  A;  i.e.,  r  =  ksrpA  where 
r ,  ksr  and  pA  are  the  reaction  rate  per  unit  area,  the  rate  coefficient 
and  the  density  of  species  A  at  the  reactive  surface,  respectively. 
They  showed  that  for  such  a  reaction,  when  a  particle  of  type  A  hits 
the  wall,  fraction  k\f  of  its  mass  is  converted  into  C,  while  the 
remaining  fraction  1  -  k\ ®  of  its  mass  is  conserved  as  A,  where 

C-(^)/(l+^)  (6) 

and  ksr  is  the  rate  constant  of  surface  reaction,  At  and  Ax  are  the  time 
and  space  intervals  in  direction  normal  to  the  surface,  and  D  is  the 
diffusion  coefficient  of  species  A.  Upon  calculating  k\f,  one  can  easily 
modify  bounce  back  boundary  condition  to  simulate  a  surface  reac¬ 
tion.  For  instance,  if  a  particle  of  type  A  moving  from  left  to  right  hits  a 
vertical  wall  on  its  path,  its  bounced  back  distribution  function  will  be 

fj  =  (i-'4B)/f 

fi  =  (i  »4B)4 

fi  =  (i-4B)/^  (7) 

/£»  =  (i-4)/A 

/?4  =  0-4B)/£ 

and  the  bounced  back  distribution  function  of  a  particle  of  type  C 
will  be 

fC  _  MWcUBfA  i  fC 

J3  —  MW/i  sr-n  1 

fC  _  MWc/.LBM  I  fC 

j  8  —  MwA/<;sr./io 

£  =  »/?+#  (8) 

fC  _  MWcUBfA  I  fC 
U3  ~  MW/i  sivll  -t-J  11 

fC  _  MW  cl, LB  f A  ,  fC 
*n4  —  MW/i  sr-n2  '  J\2 

where  MW/\  and  MWc  are  the  molar  masses  of  A  and  C,  respectively. 
Obviously  when  the  surface  reaction  is  very  fast  (/<sr  — ►  oo ),  all  the 
mass  of  a  particle  of  type  A  will  be  consumed,  thus  will 
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become  zero.  Therefore,  considering  Eq.  (7),  must  become  equal 
to  2  or  lim  k ^  =  2.  Consequently,  At  in  Eq.  (6)  must  be  chosen  as 

kST^°o 


(9) 


2.4.  Validation  of  surface  reaction  model 


Kamali  et  al.  [21]  validated  their  proposed  LB  model  by  the 
analytical  solution  of  a  one-dimensional  reaction-diffusion  prob¬ 
lem  which  is  only  valid  for  Fo  «  1  where  Fourier  number  is 
Fo  =  Dt/l2.  Since  this  condition  is  not  necessarily  valid  in  a  PEM  fuel 
cell  electrode,  we  have  validated  the  model  against  the  analytical 
solution  of  another  one-dimensional  reaction-diffusion  problem 
which  is  valid  for  all  Fo  values.  The  problem  is  a  one-dimensional 
diffusion  in  an  infinite  channel  with  two  symmetric  reactive  sur¬ 
faces  (width  =  21)  as  shown  in  Fig.  2.  It  is  obvious  that  at  the 
symmetry  line  (x  =  0),  the  density  gradient  is  zero.  Thus,  the  partial 
differential  equation  and  the  initial  and  boundary  conditions  that 
describe  the  problem  are 


pA(x,t  =  0)  =  1 

^1  =0 

8t  lx=0 

w|x=(  =  _^srP/i(X  =  1) 


(10) 


where  pa  is  the  density  of  species  A.  After  converting  Eq.  (10)  to  a 
non-dimensional  form  and  by  the  aid  of  separation  of  variables,  the 
solution  will  be  obtained  as 


P*a(x,  0 


4  sin  (An)  -X2„Fo 
2An  +  sin  (2An) 


cos 


(11) 


In  this  equation  p*A  is  the  non-dimensional  density  of  species  A , 
and  An  are  the  eigen-parameters  which  are  the  roots  of  equation 
Antan  (An)  =  Da  (given  in  Ref.  [22]),  where  Damkohler  number  is 
defined  as  Da  =  //<S/D.  The  analytical  solution  is  computed  for  the 
first  six  terms  for  several  Da  and  Fo  and  are  compared  with  the 


Fig.  2.  Scheme  of  the  problem  used  for  the  validation  of  the  modified  bounce  back 
boundary  condition. 


corresponding  lattice  Boltzmann  method  results  in  Fig.  3.  It  is 
clearly  observed  in  this  figure  that  the  lattice  Boltzmann  method 
results  are  in  a  very  good  agreement  with  the  analytical  solution  for 
all  values  of  Da  and  Fo  numbers. 


x/l 


x/l 


x/l 


Fig.  3.  Validation  of  the  modified  bounce  back  boundary  condition  by  analytical  so¬ 
lution  of  a  diffusion  problem  in  an  infinite  slab  with  reactive  surfaces  at  various  Fourier 
and  Damkohler  numbers:  (a)  Da  =  0.01,  (b)  Da  =  1,  (c)  Da  =  100;  symbols:  lattice 
Boltzmann  simulation,  solid  lines:  analytical  solution. 
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2.5.  Electrochemical  reaction 


The  electrochemical  reaction  at  the  PEM  fuel  cell  cathode  elec¬ 
trode  is  the  oxygen  reduction  reaction  as  follows 

02+4H++4e^2H20  (12) 

The  rate  of  this  reaction  per  unit  area,  r  \  is  directly  proportional 
to  the  current  density,  j,  through  r  =  j/4F  where  Fis  the  Faraday’s 
constant.  On  the  other  hand,  current  density  is  a  function  of  oxygen 
concentration  and  activation  over-potential  through  Butler— 
Volmer  equation  [2].  Therefore,  after  some  manipulation  the  rate 
equation  for  the  electrochemical  reaction  on  the  catalyst  layer, 
which  is  assumed  to  be  a  thin  interface,  can  be  expressed  as 


r"  = 


-  exp 


(13) 


where  a  is  the  roughness  factor  of  catalyst  layer,  jref  is  the  reference 
current  density,  p0,ref  is  the  reference  oxygen  density,  a/  is  the 
transfer  coefficient  for  the  forward  reaction,  ar  is  the  transfer  co¬ 
efficient  for  the  reverse  reaction,  ri  is  the  activation  over-potential 
and  p°  is  the  oxygen  density  on  the  catalyst  layer.  Eq.  (13)  can  be 
expressed  as  r  =  ksrp°  where  kSY,  the  rate  constant  of  the  electro¬ 
chemical  reaction  on  the  catalyst  layer,  is  equal  to 


(,4) 

Upon  calculation  of  /<sr,  k\ ®  can  be  calculated  for  the  electro¬ 
chemical  reaction  on  the  catalyst  layer  through  Eq.  (6)  and  then,  by 
applying  the  proposed  modified  bounce  back  boundary  condition,  the 
entire  cathode  electrode  can  be  simulated  through  an  active  approach. 


k  —  a 

Ksr  — 


f 


4 F  [pO, ref 


3.  Three-dimensional  reconstruction  of  GDL  microstructure 

In  the  early  efforts  of  lattice  Boltzmann  simulation  of  fluid  flow 
in  pore  spaces  of  a  GDL,  very  simplified  microstructures  were  used 
[23,24].  However,  with  today’s  enhancement  of  electron  micro¬ 
scopy  methods,  more  realistic  delineations  of  GDL  microstructure 
can  be  provided  by  different  three-dimensional  reconstruction 
methods.  Generally,  the  implemented  methods  for  three- 
dimensional  reconstruction  of  a  GDL  can  be  classified  into  two 
categories:  imaging  combination  and  stochastic  generation.  The 
former  captures  sequential  images  of  a  GDL  scanned  by  X-ray  or 
scanning  laser  microscopy,  and  then  integrates  these  images  to 
reconstruct  the  three-dimensional  GDL  microstructure.  The  later 
reconstructs  the  porous  medium  microstructure  based  on  the  sta¬ 
tistical  information  of  the  GDL  microstructure.  The  low  cost  and 
easy  implementation  make  the  stochastic  generation  method  a 
more  convenient  choice  for  us  than  the  imaging  combination 
method  [25].  The  GDL  microstructure  in  most  of  three-dimensional 
lattice  Boltzmann  simulations  were  reconstructed  by  stochastic 
generation  methods  [26-33].  In  the  present  study,  the  stochastic 
generation  method  proposed  by  Schladitz  et  al.  [34]  is  adopted  to 
reconstruct  three  different  carbon  paper  GDLs  which  have  also 
been  adopted  in  other  lattice  Boltzmann  simulations  [26-31].  The 
reconstruction  method  is  based  on  the  following  assumptions:  (a) 
carbon  fibers  in  GDLs  are  cylinders  having  a  fixed  and  uniform 
radius;  (b)  fibers  are  straight  and  infinitely  long;  (c)  fibers  are 
allowed  to  overlap.  The  reconstruction  method  is  started  by  sto¬ 
chastic  generation  of  a  line.  Consequently,  the  nodes  of  the  lattice 
whose  distances  from  the  line  are  less  than  or  equal  to  the  fiber 
radius  are  selected  as  solid  parts.  This  process  is  repeated  until  the 
required  porosity  is  achieved. 


The  stochastically  generated  lines  are  isotropic  in  the  material 
plane  (the  plane  of  carbon  paper)  and  are  anisotropic  normal  to  the 
material  plane.  The  degree  of  this  anisotropy  is  recognized  by  a 
factor  named  anisotropy  parameter,  /3,  which  is  0  <  /3  <  <*>  [34].  If 
0  <  /J  <  1,  the  reconstructed  carbon  fibers  are  more  probably  ori¬ 
ented  normal  to  the  material  plane,  and  if  /3  >  1,  the  reconstructed 
carbon  fibers  are  most  likely  oriented  parallel  to  the  material  plane; 
if  0  =  1,  the  reconstructed  carbon  fibers  take  a  completely  isotropic 
distribution.  For  practical  GDLs  the  anisotropy  parameter  is  a  large 
number  which  indicates  that  most  of  carbon  fibers  are  parallel  to 
the  material  plane.  Schulz  et  al.  [25]  have  reported  the  anisotropy 
parameter,  fiber  diameter  and  porosity  for  Toray090  carbon  paper 
as  10000,  7  pm  and  78%,  respectively.  In  the  present  study,  the 
Toray090  is  reconstructed  in  a  lattice  with  101  x  101  x  101  nodes 
using  the  reported  data  of  Schulz  et  al.  [25]  as  shown  in  Fig.  4(a).  In 
addition,  as  shown  in  Fig.  4(b)  and  (c),  two  other  GDLs  are  recon¬ 
structed  with  the  same  lattice  size,  fiber  diameter  and  porosity,  but 
with  different  anisotropy  parameters  of  1  and  0.01. 

4.  Computational  domain  and  boundary  conditions 

The  computational  domain  consists  of  a  portion  of  GDL,  catalyst 
layer  and  land  as  shown  in  Fig.  5.  As  mentioned  previously,  the 
catalyst  layer  is  simply  treated  as  a  thin  interface  on  the  bottom 
surface  of  the  computational  domain  on  which  electrochemical 
reaction  occurs.  Two  opposite  lateral  sides  of  computational 
domain  are  considered  as  the  inlet  (plane  x  =  0)  and  outlet  (plane 
x  =  100  pm),  while  the  other  two  lateral  sides  (planes  y  =  0  and 
y  =  100  pm)  are  assumed  as  symmetry  planes. 

The  inlet  flow  is  taken  as  dry  air  at  1.5  atm  pressure  with  oxygen 
to  nitrogen  molar  ratio  of  21/79  which  reacts  on  the  catalyst  layer 
and  produces  water  vapor.  A  constant  pressure  boundary  condition 
is  applied  on  the  outlet  such  that  a  constant  pressure  differential  of 
0.001  atm  exists  between  the  inlet  and  outlet.  To  implement  the 
constant  pressure  boundary  conditions,  the  Zou  and  He  method 
[35]  is  utilized.  Since  only  the  outlet  total  pressure  is  known  and  the 
outlet  partial  pressure  of  each  species  is  unknown,  these  partial 
pressures  are  determined  by  assuming  that  the  mole  fraction  of 
each  species  at  the  outlet  is  equal  to  that  at  the  preceding  nodes; 
this  is  a  reasonable  assumption  since  species  fraction  will  not 
change  much  through  a  lattice  unit  (1  pm).  Upon  determining 
outlet  partial  pressure  of  each  species,  constant  pressure  boundary 
condition  of  Zou  and  He  [35]  is  applied  for  each  species.  The  no-slip 
boundary  condition  is  considered  on  the  surfaces  of  all  non¬ 
reactive  solid  parts.  All  species  are  considered  as  ideal  gases  with 
different  kinematic  viscosities. 

Simulations  are  performed  on  a  desktop  computer  using  a 
parallel  FORTRAN  code  developed  in  house,  based  on  the  proposed 
model.  The  value  of  simulation  parameters  are  given  in  Table  1. 

5.  Preliminary  results  and  discussion 

The  purpose  of  this  article  is  to  present  the  details  of  the  pro¬ 
posed  three-dimensional  model  for  the  cathode  electrode  of  a  PEM 
fuel  cell.  Therefore,  only  some  preliminary  results  are  presented 
here  to  demonstrate  the  capabilities  of  the  model.  Further  simu¬ 
lations  will  be  performed  to  investigate  more  thoroughly  the  effects 
of  realistic  microstructure  of  GDLs  on  a  PEM  fuel  cell  performance. 
Such  simulations  will  require  more  sophisticated  computational 
facilities  that  are  not  currently  available  to  the  authors. 

Presented  in  the  following  are  the  simulation  results  for  the 
GDLs  described  in  Section  3  and  illustrated  in  Fig.  4(a)-(c). 

Oxygen  mole  fraction  distributions  in  the  three  reconstructed 
GDLs  are  shown  in  Fig.  6.  Cross  sections  of  carbon  fibers  are  indi¬ 
cated  in  these  figures  in  heavy  black  color.  Obviously,  for 
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Fig.  4.  Reconstructed  GDLs  with  different  anisotropy  parameters:  (a)  §  =  10000,  (b) 
0  =  1,  (c)  0  =  0.01. 

horizontally  oriented  carbon  fibers  the  cross  section  looks  more 
circular  (Fig.  6(a))  in  a  lateral  plane,  while  for  more  vertically  ori¬ 
ented  fibers  the  cross  section  looks  similar  to  an  elongated  ellipse 
(Fig.  6(c)).  It  is  also  reminded  that  the  larger  the  anisotropy 
parameter,  the  higher  probability  of  horizontal  orientation  of  the 
carbon  fibers. 

It  is  clear  in  Fig.  6  that  mole  fraction  of  oxygen  decreases  while  it 
is  penetrating  toward  the  catalyst  layer.  Oxygen  consumption  on 
the  catalyst  layer  also  causes  its  mole  fraction  to  decrease  along  the 


bulk  flow  direction  (x-axis).  However,  this  trend  is  reversed  near 
the  outlet  because  of  the  effects  of  the  GDL  pore  microstructure. 

A  comparison  among  Fig.  6(a)— (c)  shows  that  GDL  microstruc¬ 
ture  greatly  influences  species  distribution  in  the  cathode.  More 
specifically,  it  is  observed  that  when  carbon  fibers  are  mostly  ori¬ 
ented  normal  to  the  catalyst  layer  then  oxygen  density  boundary 
layer  will  be  thicker.  In  other  words,  in  this  situation  the  boundary 
layer  is  disturbed  more  strongly  by  the  fibers,  hence  a  thicker 
boundary  layer. 

Similarly,  water  vapor  mole  fraction  distributions  in  the  three 
reconstructed  GDLs  are  presented  in  Fig.  7.  Again,  due  to  its  pro¬ 
duction  by  electrochemical  reaction  on  the  catalyst  layer,  water 
vapor  mole  fraction  increases  both  normal  to  the  catalyst  layer  and 
in  the  bulk  flow  direction.  In  fact,  there  can  be  observed  a  corre¬ 
spondence  between  the  distributions  of  oxygen  and  water  vapor 
species  shown  in  Figs.  6  and  7,  since  water  production  rate  is 
directly  related  to  oxygen  concentration. 

Fig.  7(a)— (c)  illustrates  how  greatly  GDL  microstructure  can 
affect  species  distribution.  Similar  to  oxygen,  water  vapor  density 
boundary  layer  is  more  disturbed  when  carbon  fibers  are  mostly 
oriented  normal  to  the  catalyst  layer  in  a  GDL  (with  lower  anisot¬ 
ropy  parameter),  hence  a  thicker  boundary  layer.  Disturbance  of 
boundary  layer  of  water  vapor  density  may  be  regarded  as  a 
desirable  event,  since  it  can  lower  the  water  vapor  density  near  the 
catalyst  layer  which  can  mitigate  flooding  of  the  GDL.  Therefore,  the 
effects  of  GDL  pore  microstructure  on  the  form  of  water  produced 
must  be  investigated  more  thoroughly  to  help  enhance  PEMFC 
water  management. 


Table  1 

Values  of  simulation  parameters. 


Parameter 

Value 

Operating  temperature,  T 

353  K 

Operating  pressure,  P 

1.5  atm 

Pressure  differential  between  the  inlet 

0.001  atm 

and  the  outlet 

Water  vapor  mole  fraction  in  the  inlet  air 

0.0 

Oxygen  mole  fraction  in  the  inlet  air 

0.21 

Nitrogen  mole  fraction  in  the  inlet  air 

0.79 

Dynamic  viscosity  of  oxygen,  /jl° 

2.34  x  10"5  kg  nr1  s"1  [36] 

Dynamic  viscosity  of  nitrogen,  f. in 

2.01  x  10-5  kg  nr1  s-1  [36] 

Dynamic  viscosity  of  water  vapor,  /iw 

1.20  x  10  5  kg  nr1  s"1  [36] 

Oxygen  diffusivity  in  the  mixture  (Do) 

1.891  x  10"5  m2  s"1  [37] 

Reference  oxygen  concentration  (p0,ref) 

10.875  mol  m"3  [37] 

Roughness  factor  (a) 

2000  [38] 

Reference  current  density  (Jref) 

1.3874  x  10-2  A  m-2  [39] 

Transfer  coefficient  for  forward  oxygen 

0.5  [40] 

reduction  (af) 

Transfer  coefficient  for  reverse  oxygen 

1  [40] 

reduction  (ar) 
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Fig.  6.  Oxygen  mole  fraction  distribution  in  the  three  reconstructed  GDLs:  (a) 
13  =  10000,  (b)  0  =  1,  (c)  0  =  0.01. 


Fig.  7.  Water  vapor  mole  fraction  distribution  in  the  three  reconstructed  GDLs:  (a) 
0  -  10000,  (b)  0  =  1,  (c)  0  =  0.01. 


Shown  in  Fig.  8  are  current  density  distributions  on  the  catalyst 
layer  for  the  three  simulated  cases.  The  current  densities  shown 
here  correspond  to  the  “surface”  of  the  catalyst  layer.  As  mentioned 
before,  the  catalyst  layer  in  our  simulations  is  considered  as  an 
interface  between  the  GDL  and  membrane,  which  acts  as  a  charge 
conductor  itself.  The  current  density  decreases  in  the  bulk  flow 
direction  (x-axis)  due  to  oxygen  consumption  on  the  catalyst  layer. 
However,  locations  of  the  minimum  current  density  are  different  in 
the  three  cases,  which  indicate  the  influence  of  GDL  microstructure 
on  current  density  distribution. 

Although  the  variations  in  current  density  may  seem  small,  but 
one  must  realize  that  the  simulation  scale  is  also  small,  i.e.,  current 
density  variation  may  be  much  larger  for  practical  sizes  of  the 
catalyst  layer. 


An  examination  of  Fig.  8(a)— (c)  reveals  that  when  carbon  fibers 
are  mostly  oriented  normal  to  catalyst  layer  (i.e.,  smaller  anisotropy 
parameter)  a  larger  variation  in  current  density  is  observed.  This 
may  be  the  result  of  thicker  oxygen  density  boundary  layer  in  such 
a  case.  The  larger  variation  observed  in  Fig.  8(c)  matches  with  the 
lowest  overall  current  density  among  the  cases  when  carbon  fibers 
are  mostly  oriented  normal  to  catalyst  layer.  Apparently,  this  is 
caused  by  more  severe  limitation  of  oxygen  transport  to  the  surface 
of  the  catalyst  layer  due  to  the  mostly  vertical  orientation  of  the 
GDL  fibers. 

A  comparison  of  Fig.  6(a)— (c)  shows  that  the  oxygen  mole 
fraction  near  the  catalyst  layer  around  point  x  =  100  pm, 
y  =  100  pm  in  Fig.  6(a)  is  the  lowest  among  the  cases;  accordingly, 
current  density  at  this  corner  in  Fig.  8(a)  is  also  the  lowest  in 
Fig.  8(a)— (c). 
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Considering  the  results  given  in  Figs.  7  and  8,  it  is  observed  that 
when  the  orientation  of  carbon  fibers  is  dominated  in  the  hori¬ 
zontal  direction  (parallel  to  the  catalyst  layer),  a  more  uniform 
current  density  distribution  is  obtained  (desirable),  while  water 
vapor  density  is  also  greater  near  the  catalyst  layer  in  this  case 
(undesirable).  It  may  suggest  that  for  practical  purposes,  a  GDL 
microstructure  must  be  prevented  from  being  dominated  by  carbon 
fibers  oriented  in  either  horizontal  or  vertical  directions  relative  to 
the  catalyst  layer. 

Practical  GDLs  are  of  different  types  (such  as  paper,  cloth  and 
felt)  and  are  produced  by  different  companies  using  different 
techniques;  hence,  they  certainly  hold  different  microstructures. 
Moreover,  GDL  microstructure  can  be  altered  during  stack  assembly 
due  to  compression  of  the  cells.  On  the  other  hand,  the  results  from 
the  present  and  similar  studies  clearly  show  that  GDL  microstruc¬ 
ture  has  definite  effects  on  species  and  current  density  distribu¬ 
tions.  Hence,  the  microstructure  must  be  regarded  as  an  important 
characteristic  of  any  GDL  material,  and  further  research  on  its  ef¬ 
fects  must  be  pursued. 

6.  Conclusions 

In  the  present  study,  a  three-dimensional  pore-scale  model 
based  on  lattice  Boltzmann  method  is  proposed  for  the  cathode 
electrode  of  a  PEM  fuel  cell,  taking  into  account  the  electrochemical 
reaction  on  the  catalyst  layer  and  microstructure  of  the  GDL.  The 
model  enables  us  to  simulate  single-phase,  multi-species  reactive 
flow  in  a  heterogeneous,  anisotropic  GDL  through  an  active 
approach.  Reactive  flow  in  three  reconstructed  GDLs  with  different 
anisotropy  parameters  is  simulated  to  investigate  the  effects  of  GDL 
microstructure  on  species  and  current  density  distributions.  The 
results  demonstrate  that  when  carbon  fibers  are  mostly  oriented 
normal  to  the  catalyst  layer,  oxygen  and  water  vapor  density 
boundary  layers  will  be  thicker  and  more  disturbed.  Current  den¬ 
sity  variation  on  the  catalyst  layer  will  be  also  more  profound  in 
such  a  case. 

Results  of  the  present  study  display  the  important  role  of  nu¬ 
merical  methods  capable  of  taking  into  account  the  non¬ 
homogeneity  and  anisotropy  of  GDL  microstructure,  such  as  lat¬ 
tice  Boltzmann  method.  These  models  allow  one  to  hold  a  more 
clear  understanding  of  fundamental  phenomena  that  occur  in  a 
PEM  fuel  cell  electrode. 
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